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Abstract 



In this letter we apply a method recently devised in yl to find precise approximate solutions to a 
certain class of nonlinear differential equations. The analysis carried out in 1] is refined and results 
of much higher precision are obtained for the problems previously considered (Duffing equation, 
sextic oscillator). Fast convergence to the exact results is observed both for the frequency and 
for the Fourier coefficients. The method is also applied with success to more general polynomial 
potentials (the octic oscillator) and to the Van Der Pol equation. 
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I. INTRODUCTION 

In this letter we extend the analysis of Q], where a new method for the solution of os- 
cillatory nonlinear problems was devised. This method works by combining the Lindstedt- 
Poincare (LP) method to the Linear Delta Expansion (LDE) , or Optimized Perturba- 
tion Theory (OPT) 4]. A detailed list of references on the LDE can be found in 

The LP method is used to solve approximately nonlinear differential equations with pe- 
riodic solutions, by introducing strained coordinates, which allow to eliminate the secular 
terms arising in the perturbative expansion; in the LDE method an insoluble theory is in- 
terpolated with a soluble one depending on some arbitrary parameter: an expansion is then 
carried out in terms of a suitably defined "perturbation". Since the "perturbation" is not 
a priori fixed and not necessarily expressed in terms of "small" parameters, the method is 
truly nonperturbative. 

As proved in the application of LDE to the LP method, which is at the core of our 
method, allows to extend the analysis to the nonperturbative regime (large nonlinearities) 
and to obtain errors which are small even in regimes where the LP method fails completely. 

In this letter we pursue two goals: on one hand to investigate the convergence of our 
method in the cases previously analyzed to higher orders; on the other hand to test the 
method by applying it to more demanding problems, such as the Van der Pol equation. 

The letter is organized as follows: in Section |Tll we study the Duffing equation and discuss 
the convergence of the method; in Section lTTTl we analyze more general anharmonic potentials 
(sextic and octic); finally, in Section HVI we consider the Van der Pol equation. Finally in 
Section El we draw our conclusions. 

II. DUFFING EQUATION 

We start our analysis by considering the Duffing equation, which describes the oscillations 

2 .4 

of a unit mass in a potential V{x) = ^ + ^^-j-- It reads: 

^{t) + x{t) = -^^x^t). (2.1) 
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/i is a coupling which controls the strength of the nonlinear term. The application of the 
LPLDE method requires to write eq. ()2.H) as 

^ + (l + = ^ [-^ ^'(^) + ' (2-2) 

where A is an arbitrary parameter and Q is the exact frequency of the system. We also 
introduced a power-counting parameter 6, which allows to recover eq. ()2.1|) when the value 
5 = 1 is taken. The r.h.s. of eq. ()2.2|) is then treated as a perturbation, although its size is 
not fixed, being A completely arbitrary. When the perturbative expansion is carried out to 
finite order, a spurious dependence upon this parameter will show up in all the observables. 
The Principle of Minimal Sensitivity (PMS) will then be used to choose a A which minimizes 
this dependence. This procedure is explained in detail in j]]. 

To a given order Nmax we write the frequency and the solution as 

n'= , x{t) = E '^^ ^n(r) . (2.3) 

n=0 n=0 

and we take 5 = 1 aX the end. Both the coefficients a„ (which are obtained by eliminating 
the resonant contributions at any given perturbative order) and the solutions a;„(r) are A- 
dependent. The dependence upon the arbitrary parameter is then eliminated by imposing 
the PMS condition = 0. Remarkably we have observed that the solution to this equation 
to any given order is approximated extremely well by the third order solution, i.e. A = 
^ ^ both for positive and negative values of /i. This allows us to obtain expressions 
which are completely analytical. We can also attach a more physical meaning to the PMS 
condition by looking at Fig. ^ in this Figure we plot the energy of the oscillator evaluated 
at x = 0, where all the energy is kinetic, and compare it with the exact expression for 
the energy, which is given by the potential energy at the inversion point. The approximate 
solution found with eq. ()2.3p violates the conservation of energy by an amount which depends 
upon A. However, when the optimal value of A is used (the vertical line in the Figure) such 
violation is nearly minimal. 

We have noticed that the coefficients a of eq. ()2.3j) can be written as: 

3 .n l^2n (^^ ^A^^ 

aQ=l + - A ^ , a2n= , tt2n+l = (2.4) 

4 (l + f ^2 ^) 

where the k„ are purely numerical coefficients and n = 1,2,.... Such relations hold for 
all the a that we have calculated. In Fig. |21 we plot the logarithm of the coefficients, for 
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FIG. 1: Energy of the Duffing oscillator corresponding to x = O.The vertical line is the optimal 
value to third order. We assume /i = 100 and A = 1. 
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50; the dashed line is a linear fit corresponding to: 



Hn = 0.0663 ■ e 



-1.46225 n 



(2.5) 



with n even. The first few coefficients are written in Table HI 
The expression in eq. f|2.3|) can be written equivalently as 



^ max 



X 



approx 



it) = E 



c. 



(approx) 



COS [{2n +1) nt] . 



(2.6) 



n=0 



where qI^pp'^"^) are the approximate Fourier coefficients obtained with our method. As for 
the an, also the c^^pp"^"^) are analytical and take the form 

ax 



^(approx) 



m=0 



(2.7) 



where Cnm are the corrections of order m to the Fourier coefficient corresponding to 
cos [(2n + 1) Q t]. The latter can be written as: 

- _ Pnm ^ (^^ f^) /r, o\ 

Cnm — / q ,n \" ' 
(l + I /i) 

where /Snm is a numerical coefficient. In Fig. El we display the logarithm of /3 corresponding 
to different Fourier modes and to different orders in our method: interestingly, we observe 
that also the f3s decay exponentially with the order of the expansion. 
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FIG. 2: Coefficients of the formula H2.5|l . The dotted Une is a Unear fit. 



TABLE I: First few coefficients of the expansion of eq. (|2.4|) . 



n 




n 




2 


3 
128 


12 


19974549 
18014398509481984 


4 


51 


14 


128255751 


131072 


2305843009213693952 


6 


213 


16 


435036452211 


16777216 


151115727451828646838272 


8 


70515 


18 


2950668677535 




137438953472 


19342813113834066795298816 


10 


406179 


20 


163068192461619 


17592186044416 


19807040628566084398385987584 




FIG. 3: Logarithm of /? corresponding to different Fourier coefficients and to different orders. 
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We now present some results obtained with our method. We will first consider the case 
of positive //. In Fig. |3]we plot the logarithm of the error defined the equation 

X 100 , (2.9) 

where fl^ is the squared frequency obtained with the LPLDE method and fll^act is the exact 
value. The blow-up in the figure allows to better appreciate the small differences in the 
error corresponding to the different values of /i (we use /i = 10, 100, 10000). We notice that 
the error is practically unaffected by the size of /i. This result is clearly understood by eq. 
()2.4|) . because the coefficients a„ clearly go to faster than the for /i > 0. 

In Fig. El we compare the coefficients c^^pp*""^) of the approximate solution with the ones 
of the exact (numerical) solution: 

oo 

Xexactit) = Ci^""*^ COS [(2n + 1) ^] t] . (2.10) 

n=0 

Notice that the approximate series ()2.6|) is truncated at the maximum frequency {2Nmax + 
1) ^mo.x] however this cutoff frequency can be increased by applying the method to higher 
orders |7|. 

Clearly our method is capable of reproducing the first few coefficients of the Fourier 
series with great accuracy. Although the modes with higher frequency turn out to be poorly 
approximated, they don't affect the overall quality of the approximation, given the small 
size of their contributions. 

We now consider fi = —1. In this case the potential has maxima located at a; = ±1 and 
an oscillatory behaviour is permitted only for amplitudes A < 1. A = ±1 are points of 
(unstable) equilibrium and therefore the period diverges in correspondence of these values. 

By looking at Fig. IHl we see that the slope of the logarithm of the error is now strongly 
dependent upon the amplitude, in contrast with the case previously analyzed. This result 
can be understood in view of eq. ()2.4|) : because A"^ // is now negative, the sign in the equation 
is changing order by order and the size of the denominator is smaller. 

In Fig.[7|and|Hlwe analyze the Fourier coefficients of the approximate and exact solutions 
and then plot the difference between the two as a function of time. We consider oscillations 
with an amplitude very close to the maximum, i.e. A = 0.99. Also in this case we notice 
that the approximation works very well. 
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FIG. 4: Logarithm of the percentile error defined in eq. (|2.9|) for A = 10 and different values of 
(/" > 0). 
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FIG. 5: Left plot: exact (numerical) and approximate coefficients of the Fourier series corresponding 
to A = 10 and fi = 10^; Right plot: ratio of the approximate to the exact coefficients. 



III. ANHARMONIC POTENTIALS 



The analysis carried out in the previous Section can be extended easily to more general 
anharmonic potentials of the form: 

V{x) = -x' + -^x^'' (3.1) 

with positive integer. In particular we study here the sextic and octic oscillators: we will 
see that many of the features observed previously for the Duffing equation hold true also in 
this case. Remarkably the optimal value of A calculated to third order falls very close to the 
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FIG. 6: Logarithm of the percentile error defined in eq. (|2.9j) for /u = — 1 and different amplitudes. 
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FIG. 7: Left plot: exact (numerical) and approximate coefficients of the Fourier series corresponding 
to A = 0.99 and /.i = —1; Right plot: ratio of the approximate to the exact coefficients. 

exact solution of the PMS condition to any given order. This result will allow us to obtain 
fully analytical expressions also in the case of sextic and octic oscillators. 
We obtain indeed the expressions for the optimal values of A: 
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By looking at Fig. Inland ^1 we see that our method provides an excellent approximation 
both in the regime of positive and negatives /i: indeed we see that the error decays expo- 
nentially following a law which is very similar to the one observed in the Duffing oscillator. 
These similarities reflect also in the behaviour of the Fourier coefficients, which however we 
will not display here. 
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Left plot: exact (numerical) and approximate solutions of the Duffing equation corre- 
',to A = 0.99 and ^ = — 1; Right plot: the difference Xexactit) — Xapprox{t)- 
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FIG. 9: Left plot: logarithm of the error over the squared frequency for the sextic oscillator with 
/i = 1, 100, 10000. Right plot: logarithm of the error over the squared frequency for the sextic 
oscillator with A = 0.5,0.9 and 0.99 and /u = —1. 

IV. VAN DER POL EQUATION 

We now come to consider the Van der Pol equation (see for example P]): 



(4.1) 



This equation possesses a limit cycle and can sustain periodic oscillations with definite 
period and amplitude depending upon the constant However, unlike the cases treated 
before, eq. ()4.1|) does not correspond to a conservative system and indeed the r.h.s of the 
equation either damps or enhances the oscillations depending upon the size of x. 
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FIG. 10: Left plot: logarithm of the error over the squared frequency for the octic oscillator with 
/i = 1, 100, 10000. Right plot: logarithm of the error over the squared frequency for the octic 
oscillator with A = 0.5,0.9 and 0.99 and /u = —1. 
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We will now apply our method to this problem in the usual fashion by writing: 
x+(l + X^) x{t) = 6[hQ (l-x^) x + X^ 
where the expansion 

= E 7n , X{t) = E ^n{r) 



(4.2) 



(4.3) 



n=0 n=0 

is assumed. The coefficients 7„ are fixed by the requirement that the resonant terms at a 
given order n vanish. 

As we can see from the left plot of Fig. ^2 the optimal value of A turns out to be strongly 
order dependent for a given /i unlike in the cases considered previously. However, at a fixed 
order, we observe that the optimal value of A grows linearly with ^ (right plot). The dashed 
line in the right plot is the fit A = 0.212599 + 1.17166 This behaviour of A allows to 
obtain a semi-analytical formula both for the period and for the solution which works for 
moderate values of /i, where the LP method is already not applicable. 

In Table ITTl we report the period of the solution of the Van der Pol equation as obtained 
in our method (to order 44) for /i ranging from 1 to 10 and we compare it with the numerical 
results of Notice that, at the largest value of /i considered here (/x = 10), we obtain an 
error of 12 %. Much smaller errors are obtained at lower values of /i. 

In Fig. ^1 we display the approximate (to order 43) and exact (numerical) solutions of 
the Van der Pol equation for = 3. For such a large value of // our method provides an 
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FIG. 11: Left plot: Optimal value of A as a function of the order for /i = 3. Right plot: Optimal 
value of A to order 44 as a function of /x. The dashed line is the fit A = 0.212599 + 1.17166 /i. 



TABLE II: Comparison between the approximate and exact period of the Van der Pol equation. 





T 

approx 


'^exact 




T 




1 


6.66328685 


6.66328686 


6 


12.83944539 


13.06187474 


2 


7.62995604 


7.62987448 


7 


14.03245644 


14.53974774 


3 


8.86085271 


8.85909550 


8 


15.09907603 


16.03817623 


4 


10.19701347 


10.20352369 


9 


16.02685408 


17.55218414 


5 


11.54602911 


11.61223067 


10 


16.8128186 


19.07836957 



excellent approximation both to the period and to the function. Notice that in this regime 
the LP method is not applicable. 



V. CONCLUSIONS 

In this work we have extended the analysis carried out in Q and we have applied the 
LPLDE method to much higher orders: we have observed that the results obtained converge 
quickly to the exact values, with errors which decay exponentially with the perturbative or- 
der. This behaviour was observed for the Duffing equation, for the sextic and octic oscillator. 
Remarkably the results obtained to any given oidei with our method are fully analytical. 

The method was then tested on the Van der Pol equation where the calculations were 
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FIG. 12: Approximate and exact solutions of the Van der Pol equation to order 43 for = 3. 
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FIG. 13: Approximate and exact Fourier coefficients of the Van der Pol equation calculated at 
order 43 for fi = 3. Left plot: coefficients of the cos; Right plot: coefficients of the sin. 

performed up to order 44. In this case, we were able to obtain good results for moderate 
values of f« 3), where the LP method already fails. These results suggest that the 
extension of the method to higher orders could extend the region of application to larger 
values of fi. Moreover, although the optimal value of the variational parameter is not known 
analytically, as in the previous case, the use of a /x dependent fit allowed us to obtain also 
in this case a semi-analytical approximation for the frequency (period) which holds in the 
region of moderate /i. 

It remains to explore the application of the present method to a wider class of nonlinear 
problems with periodic solutions, as well as studying the Van der Pol equation to higher 
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orders. This will be done in future work. 

VI. ACKNOWLEDGMENTS 

The Authors would like to thank A.Aranda and R. Saenz for useful discussions. They 
also acknowledge the support of the "Fondo Alvarez- Buylla" of the University of Colima 
and of Conacyt grant no. C01-40633/A-1. 



13 



[1] P.Amore and A.Aranda, Phys. Lett. A 316, 218-225 (2003) 

[2] A. Lindstedt, Mem. de I'Ac. Imper. de St. Petersburg 31, 1883 

[3] A. Okopmska, Phys. Rev. D 35, 1835 (1987); A. Duncan and M. Moshe, Phys. Lett. B 215, 
352 (1988) 

[4] V. I. Yukalov, Phys. Rev. A 58, 96 (1998); J. Math. Phys. 32, 1235 (1991); Teor. Mat. Fiz. 28 
(1976) 92. 

[5] A. Buonomo, SIAM Journal of AppUed Mathematics, vol. 59, No. 1, pp. 156-171 

[6] M. Strasberg, Recherche de solutions periodiques d'equations differentieUes non hneaires par 

de methodes de discretisation, P. Jannsenns, J. Mawhin and N. Rouche, ed. Hermann, Paris, 

1973, pp. 291-321 
[7] The results displayed in the Figure are calculated to order 50. 



14 



